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Abstract. Many real-world complex systems are best modeled by multiplex networks of interacting network 
layers. The multiplex network study is one of the newest and hottest themes in the statistical physics of 
complex networks. Pioneering studies have proven that the multiplexity has broad impact on the system’s 
structure and function. In this Colloquium paper, we present an organized review of the growing body of 
current literature on multiplex networks by categorizing existing studies broadly according to the type of 
layer coupling in the problem. Major recent advances in the field are surveyed and some outstanding open 
challenges and future perspectives will be proposed. 

PACS. 89.75.He Networks and genealogical trees - 89.75.-k Complex systems 


1 Introduction 

Many real-world complex systems ranging from living or¬ 
ganisms and human society to transportation system and 
critical infrastructure operate through multiple layers of 
distinct interactions among constituents as well as the 
interplay between these interaction layers to fulfill their 
emergent function [T1I21I3] . Multiplex network [3Tl5ll51f7U 51 
I91IT01HT] is a class of networks introduced to better model 
such systems, in which the same set of nodes are connected 
via more than one type of link^] (Fig. 1). Each type of 
links in multiplex networks constitutes the network layer. 
Examples for multiplex networks abound: Individuals in 
a society are networking through numerous social rela¬ 
tionships such as friendship, kinship, co-workership and 
via a multitude of communication channels such as online 
and offline contacts [12]. Critical infrastructure provides 
essential support for the functioning of modern society 
through concerted operations of multiple interlinked and 
interrelated networks such as energy production and sup¬ 
ply, telecommunication, and transportation networks H3|. 
The study of multiplex networks has emerged as one of 
the major contemporary topics of network theory HUSH 
Eumumim, along with the parallel development of closely 
related topics such as interacting networks [14], interde¬ 
pendent networks [15], and interconnected networks [RE 
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1 In the literature, the term “multiplex” tends to be used in 
a more loose sense, referring also to more general multilayered 
and interconnected systems. Throughout this paper, however, 
we will try to adhere to the strict definition of the term. 



Fig. 1 . An illustrative example of the multiplex network 
of nine nodes with two layers, the red (solid) and the blue 
(dashed) layer. A node and its replica are shown linked by the 
dotted line, to denote the identity relation. In addition to the 
multiplex network topology like this, the type of layer coupling 
needs also to be implemented correctly, for a reliable treatment 
of given multiplex problem. 


Introduction of manifold layers in multiplex networks 
necessitates new conceptual as well as computational de¬ 
velopments beyond the standard single-network theoretic 
framework that has been well-established during the last 
fifteen or so years m- The multiplex network is not just a 
simple stack of many network layers. What makes a truly 
multiplex system is the functional coupling between layers, 
which gives rise to in general non-additive and nonlinear 
effect for the emergent function of the system. Hence a 
proper consideration of the contexts of multiplexity that 
the multiple layers make is essential for correct description 
and successful understanding of multiplex problems. 
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There are two important facets of multiplexity, which 
we might refer to as (i) the pattern of multiplexity and 
(ii) the type of multiplexity. The former concerns how 
the layers are coupled structurally. The same set of net¬ 
work layers can be coupled in many ways to form different 
multiplex structures (and vice versa). In real-world mul¬ 
tiplex systems, the coupling structure between layers are 
far from random but correlated, a property termed the 
correlated multiplexity [5], which has been assessed to be 
significant empirically [18] in terms of the interlayer de¬ 
gree correlation and the link overlap. The type of multi¬ 
plexity concerns how the layers are coupled functionally , 
that is, how the function of one layer affects that of an¬ 
other (and vice versa). For instance, the multiple layers 
in a multiplex system may be coupled either in comple¬ 
mentary/connective manner as in the transportation sys¬ 
tem jHimni] or in cooperative/dependent way as in the 
critical infrastructure dnni5]. Depending on the type of 
multiplexity, same multiplex structures can behave quite 
differently. A legitimate understanding of multiplex net¬ 
work problems can thus be achieved only when both struc¬ 
tural and functional aspects of multiplexity are correctly 
taken into account. 

In this Colloquium paper, we attempt an organic 
overview of the fast-growing body of current literature on 
statistical physics of multiplex networks. We will pay par¬ 
ticular attention to put the diverse approaches in related 
problems into a coherent setting characterized and catego¬ 
rized broadly by the type of multiplexity. In so doing, we 
hope this paper not to reiterate but to complement several 
existing review articles on closely-related topics. Readers 
are referred to references mm for more focused sum¬ 
mary on interdependent networks, to references [2,3, for 
broader expositions on multilayer networks, and to refer¬ 
ences mm for comprehensive reviews on temporal net¬ 
works. A compendium volume [1] collecting contributions 
from pioneering groups is also available, providing an over¬ 
arching overview of the early development of the field at 
large. 

The rest of this Colloquium paper is organized as fol¬ 
lows. In Section 2, we discuss how to characterize and 
model the multiplex systems. There will also be intro¬ 
duced the basic definitions and terminology. In Section 3, 
we review multiplex problems with cooperatively-coupled 
layers, such as mutual percolation and multiplex cascades, 
for which major new concepts and analytic methods spe¬ 
cific to multiplex systems are introduced. In Section 4, 
problems with complementary layers, such as transport 
and epidemic spreading on interconnected layers, are sur¬ 
veyed. In Section 5, other types of layer coupling such as 
competitive coupling and directional coupling are briefly 
discussed. The paper will close with conclusion and out¬ 
look in Section 6. 

2 Characterization and modeling of the 
multiplex structure 

A first step towards understanding of new network con¬ 
cept is to formulate appropriate means to characterize the 


structural properties. To motivate, we first look for em¬ 
pirical data for multiplex networks, ranging from social 
and man-made systems to biological systems. Basic rep¬ 
resentations of multiplex network ensemble and various 
structural measures are introduced. Finally, random graph 
models of multiplex networks will be introduced as theo¬ 
retical tools for systematic and analytic understanding of 
the multiplex system. 


2.1 Empirical studies and multiplex network data 

Empirical studies on real-world network data have been 
a major fuel of network science, motivating newer and 
more sophisticated network concepts and developments of 
necessary analytic tools. 

Society is a prime example of the multiplex network 
system. Since early social network studies, data on the so¬ 
cial multiplex (sometimes also termed as multi-stranded) 
relationships had been collected through surveys and in¬ 
terviews, and its social implications had been discussed [251 
l24U25i[2(T0 With the recent development of digital tech¬ 
nology, human individual activities are being recorded in 
unprecedentedly high spatiotemporal resolution and in large 
scale. Studies exploiting this opportunity came from the 
analyses of the massive multiplayer online games nmzi 
[28] . first quantifying the correlated multiplexity in large 
scale. Another kind of large-scale yet more accessible mul¬ 
tiplex social network data is the communication-collaboration 
network between scientists [2911301131] . 

Many man-made systems in modern society are also 
multiplex network systems. An example is the transporta¬ 
tion system. In the scale of urban and nation-wide tran¬ 
sit system, several datasets on different transportation 
modes [321133] have been constructed. At a larger scale, 
a multiplex European airway network with thirty-seven 
layers corresponding to routes operated by different com¬ 
mercial airline companies [34] and the worldwide airport- 
seaport coupled network [35 have been analyzed. An¬ 
other example is the so-called critical infrastructure sys¬ 
tems ranging from the power grid and water supply net¬ 
works to wired and wireless information communication 
infrastructures HEZ], which motivated the concept of 
interdependency between layers in a multiplex system, 
leading to interesting new physics (see Sect. 3 for more 
details). The global economic system is also a multiplex 
system consisting of various financial and political layers. 
Best characterized example is the international trade net¬ 
work among countries as a multiplex network with layers 
representing different types of goods [2139]. 

Biology is governed by multiplex networks across many 
levels of hierarchy. At the cellular level, the basic cellular 
function is fulfilled by the coordinated actions of many 
layers of biomolecular interactions. Reconstruction and in¬ 
terrogation of major biomolecular network layers in model 
organisms gnum and human [431144] have made steady 

2 In the social network literature the term multiplexity has 
often been used to refer the degree of overlap in relations (links) 
across different layers. 










Kyu-Min Lee, Byungjoon Min, and K.-I. Goh: Towards real-world complexity: an introduction to multiplex networks 3 


progress in the last decade towards the whole-cell model¬ 
ing [45] . At the physiological level, an example is provided 
by the neuronal network of Caenorhabditis elegans consist¬ 
ing of two layers of connections, the synaptic connections 
and the gap junctions [46], and the layer-level structural 
properties have been analyzed in reference [47]. At the 
ecological level, species interact with one another through 
mutualistic, host-parasite or predator-prey relationship, 
which should be treated in terms of the multiplex net¬ 
work [48] . 


2.2 Mathematical formulations and measures 

On first thought, it might seem straightforward to general¬ 
ize various network concepts and measures well-defined for 
single-layer networks PH into multiplex networks by intro¬ 
ducing additional index for the network layer. This turned 
out not always to be the case [49]. Often such simple- 
minded generalization could miss the essential feature of 
the multiplex system, that is the context of multiplexity. 
In this section, we shall try to summarize the current sta¬ 
tus of theoretical development in this direction, clarifying 
what has been done and what has to be done. For the 
sake of simplicity of presentation, we shall consider only 
the multiplex networks consisting of layers of simple undi¬ 
rected unweighted graphs. 

Throughout this paper, we will use the term simple 
multiplex to refer the multiplex networks in which a node 
participating to multiple layers (to be called a multiplex 
node) should participate to all the layers. This restriction 
is relaxed in more general multilayer interconnected net¬ 
works. When every node in the network is multiplexed, 
the system is said to be fully multiplexed. Otherwise, it is 
called partially multiplex. 


2.2.1 Matrix representation 

One can fully represent the given instance of network 
(or graph) with the adjacency matrix A with elements 
Aij = 1 if the nodes i and j are connected by a link and 
Aij = 0 otherwise. A straightforward generalization of the 
adjacency matrix for multiplex network is the so-called 
supra-adjacency matrix m or adjacency tensor [4]. For a 
multiplex network with N nodes and £ layers, the supra- 
adjacency matrix A is an £ x £ matrix of N x N blocks 
of layer-to-layer adjacency matrices. The (ij)-element of 
(afr)-block, (A a b)iji fo r a simple multiplex is thus given by 

, . x _ J 1 , if nodes i and j are linked within layer a, 

{ aajij | o, otherwise, 

(1) 

and 

, , x _ J 1 , if node i is present in both layers a and 6, 

{ ab)m | q, otherwise, 

( 2 ) 


with i, j = 1, • • • , AT, and a, b = 1, • • • , A Throughout this 
paper, the alphabetical indices (a, 6,...) will be used for 
the dummy index of layers, with the indices i and j re¬ 
served for the node index. In the strict definition of sim¬ 
ple multiplex networks, the off-diagonal blocks are always 
diagonal matrices. Therefore the supra-adjacency matrix 
representation is somewhat redundant to encode the link 
structure of multiplex networks, but we shall reserve this 
representation for its applicability to more general mul¬ 
tilayer interconnected networks in which the interlayer 
links can connect different nodes across layers, resulting in 
nonzero off-diagonal elements in the off-diagonal blocks ISO! 
Initial idea had been applied to study the community 
structure [4]. Supra-Laplacian matrix of the multiplex net¬ 
work can also be defined by following the similar gen¬ 
eralization procedure [7], which has been applied in for¬ 
malizing diffusion-type linear dynamic processes on mul¬ 
tiplex networks [7]I51] I52 ] (see Sect. 4.2 for more details) 
and in generalizing the eigenvector centrality |53|. Gen¬ 
eralization of the celebrated PageRank algorithm to the 
multiplex network has also been proposed [54]. Despite its 
mathematical elegance and theoretical appeal, the supra- 
adjacency matrix based representation per se is not suffi¬ 
cient in that it does not fully encode the functional context 
of multiplexity, which should be specified separately. 

2.2.2 Multiplex degree 

The multiplex degree of node i (or multidegree as is some¬ 
times called) is constructed as a list of its degrees in each 
layer, 

ki = {*i 1) ,... ) ^} 1 (3) 

with k \ a ^ = J2j(A aa )ij, where the layer index is writ¬ 
ten superscripted to avoid confusion with the node index, 
(note that in the general multilayer case, interlayer degrees 

{k \ ah )} need also be specified). 

The concept of degree centrality can be generalized in 
many ways based on the multiplex degree. Proper defi¬ 
nition of degree centrality would depend on the system’s 
functional characteristics and the context of multiplex¬ 
ity. For example, in a system with complementary layers 
the total degree of a node could be a good candidate for 
a degree centrality measure. In the presence of overlap¬ 
ping links, however, either the total sum of the degree 
in each layer or the total number of non-redundant links 
over all the layers could be more meaningful, depending 
on the physical process for which the centrality aims to 
address. Moreover, in systems with cooperative or com¬ 
petitive layers it would not even be clear whether the total 
degree could be a meaningful quantity for centrality at all 
(considering, for example, the case of mutual connectivity, 
which is discussed in detail in Sect. 3). In general, the con¬ 
cept of centrality of a node should be generalized with care 
in multiplex networks, since a context-blind generalization 
can lead to irrelevant or even misleading implications. 

The joint distribution of multiplex degree P(k) is a 
natural generalization of the degree distribution P(k) for 
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Fig. 2. Schematic illustrations of three patterns of interlayer degree-correlated multiplex networks proposed in 0, the maximally 
positive (MP), uncorrelated, and maximally negative (MN) cases of two layers (the red solid and the blue dotted layers). 


single-layer network, encoding fundamental structural in¬ 
formation of a multiplex network ensemble D3I- 

2.2.3 Correlations between layers 

In multiplex networks, the correlation property of the layer 
coupling is an important part of the context of multiplex- 
ity, introducing extra dimension of network correlations in 
addition to the correlation properties of individual layers. 
The lowest-order layer-level correlation would be the ex¬ 
tent of node multiplexity, viz., what fraction of nodes in 
the network appear across different layers (multiplexed). 

To the next order, the interlayer degree correlation 
specify how correlated the degrees of a node are across dif¬ 
ferent layers. It is easy to conceive the prevalence of this 
type of correlation in multiplex systems: a person with 
many friends would likely have many colleagues also in 
the workplace, being a friendly person; an airway-hub city 
tends to be very well-connected in ground as well, through 
railways and highways. In the presence of interlayer de¬ 
gree correlation the joint degree distribution P(k) does 
not factorize into the product of individual layer’s degree 
distribution. In parallel with the assortativity measure for 
single-layer network [55], simplified measures based on the 
correlation coefficients were introduced to quantify the in¬ 
terlayer degree correlation: Pearson correlation coefficient 
between two layers a and 6 , 

(k a kb) {ka)(kb) / A \ 

Pab — 5 (4) 

a k a CTk b 

where k a ^ denote the degrees of the same nodes in differ¬ 
ent layers, (x) the average over all multiplex nodes, and 
a x the standard deviation thereof, has been used for this 
purpose in empirical [18] as well as theoretical studies [9]. 
Other correlation coefficients like Spearman rank correla¬ 
tion or Kendall’s r rank correlation have also been consid¬ 
ered m- The average degree in layer a, (k a ), of nodes hav¬ 
ing the same degree kb in the other layer b , (fc a )(fc&), can 
also encode the interlayer degree correlation property: In¬ 
creasing (decreasing) behavior of ( k a )(kb ) signifies positive 


(negative) interlayer degree correlation [56]. More kinds 
of interlayer correlations were introduced and empirically 
measured in reference m • 

To investigate systematically the effect of interlayer 
degree correlation, it has been proposed to compare three 
specific patterns of correlated coupling, the maximally- 
positive (MP) coupling, random (uncorrelated) coupling, 
and the maximally-negative (MN) coupling |5] (Fig. 2). 
Given the two layers, the MP-coupled multiplex is con¬ 
structed by matching the nodes of the same degree-rank 
from each layer. That is, the largest (lowest) degree node 
in one layer is coupled to the largest (lowest) degree node 
in the other layer. The MN coupling is constructed in the 
opposite way of MP coupling. Effects of interlayer de¬ 
gree correlation have been studied in this way for per¬ 
colation [5], mutual percolation and robustness 153. and 
cascading failures [58] problems. 

Another important layer-level multiplex correlation is 
the presence and extent of link overlap across different 
layers. In social network literature, the term multiplexity 
often refers to this feature, asking questions like how the 
existence of link in one layer facilitates or constrains the 
formation of the link between the same node pair in an¬ 
other layer and what the distinctive functional role of the 
overlapping link is to various social dynamical processes. 
From theoretical point of view, the link overlap would be 
statistically unlikely to exist if the layers were coupled 
completely randomly (assuming sparseness of the layers). 
Therefore the presence of link overlap signifies underlying 
non-randomness of layer coupling in the system. Empiri¬ 
cal studies have shown that the link overlap between two 
social network layers is indeed very significant, with Jac- 
card index as high as 0.16 in the example of Pardus net¬ 
work [18 . Functional roles of overlapping links have been 
studied theoretically for mutual percolation and 

viability m , revealing the intricate role that overlapping 
links play in cooperatively-coupled multiplex networks. 

To a higher-order correlation, the closed three-body 
relation, or the triad, can measure the degree of transi¬ 
tivity of interactions, and has served as the basis for the 
notion of network clustering m As triangular relations 
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in multiplex networks can be closed not just within a layer 
but also across layers, generalized notion is required. For 
example, mutual friends of the same person can form busi¬ 
ness relationship more easily than other random pairs, 
through the brokering action by the mutual friend. The 
so-called cross-layered clustering coefficient has been de¬ 
fined (621163i[6ffil65j to address this new type of transitivity 
in multiplex networks. A family of clustering coefficients 
parametrized with the layer coupling strength has also 
been proposed [66! . The functional meaning and impact of 
cross-layer transitivity for cooperatively or competitively 
coupled layers (such as mutual percolation introduced in 
Sect. 3) remain to be fully examined. 

2.2.4 Shortest path and other distance-based measures 

Similar conceptual complication also applies to the no¬ 
tion of shortest paths in multiplex networks and related 
quantities like pairwise distance and centrality measures 
like closeness, betweenness, and efficiency PH- For sys¬ 
tems with complementary/connective coupling (imagine 
for example the interconnected transportation network) it 
is relatively easy to generalize the notion in terms of the 
optimal paths along which one can reach from one node 
to another in minimum time, by incorporating the differ¬ 
ence in link capacity (say different average speed of each 
transportation layer) and the layer-switching cost (say av¬ 
erage transit time). Quantities defined along this line have 
been used in transport problems on multiplex networks 
|62ll65il67j . Still, the meaning and relevance of the shortest 
path and distance-based measures in systems with cooper¬ 
atively or competitively coupled layers are much less clear 
(as in the case of mutual connectivity in Sect. 3). Proper 
distance measures on multiplex networks should thus be 
context-dependent and incorporate the functional aspect 
of the layer coupling and the dynamic process under con¬ 
sideration, calling for more attention. 

2.3 Multiplex network models 

2.3.1 Multiplex random graph models 

Random graph models play two important roles in net¬ 
work theory. First they provide null models for network 
structure with which empirical structure of real-world net¬ 
works can be compared and key structural characteristics 
can be extracted. Second they serve as theoretical testbeds 
on which systematic analytic and computational analyses 
can be performed to establish basic understanding of the 
structural and dynamical characteristics of the networks. 

Random graph ensembles are constructed by specify¬ 
ing structural constraints, such as the number of nodes 
and links, degrees of nodes, degree-degree correlations be¬ 
tween connected nodes, etc ., as one tries to model net¬ 
works with higher order correlations. Most elementary and 
widely-used graph ensemble is the one with the degree se¬ 
quence or degree distribution constrained. Such an ensem¬ 
ble is straightforward to generalize to multiplex networks 


by using the multiplex degree k = {k a } or its distribution 
P(k) [5l l8lll4[|15ll29] . (In more general multilayer case, the 
interlayer degrees {/q^} need also be used |8lfl4] b Using 
such random graph models, percolation problems of un¬ 
correlated interconnected M and interdependent layers 
m were first studied. Percolation of multiplex networks 
with interlayer degree correlation was investigated in refer¬ 
ences j5ll29l[35ll571l68] . It is worthwhile to note also the re¬ 
lated predecessors like random graphs with colored edges 
[69ll7Qj and the multi-component static scale-free network 
model [71] . 

Networks with finite fraction of overlapping links can¬ 
not be constructed in the above way. To handle this, one 
need to constrain the degrees for the overlapping links sep¬ 
arately. In the case of multiplex networks with two layers, 
for example, this can be done by decomposing the mul- 
tiplex degree as k = {k A ,k B ,k overlap}, with Overlap rep- 
resenting the degree for overlapping links and the other 
two for non-overlapping links. The random multiplex en¬ 
semble with link overlap was introduced in 0 (th erein 
the overlapping links were called the multilinks). Using 
this method the mutual percolation property of multi¬ 
plexes with link overlap were investigated [5911601161] ■ Spa¬ 
tially embedded multiplex ensemble was also proposed as 
a possible way to produce the link overlap in multiplex 
networks m- 


2.3.2 Growing multiplex network models 

Growing network models represent another important class 
of network models, in which the generative rules specifying 
how the new link is formed as the new nodes and links are 
introduced in (pseudo) time to the network define the par¬ 
ticular random graph process. This class of models have 
been useful in providing insights on the relationship be¬ 
tween the microscopic linkage rules and the macroscopic 
network properties m, as in the famed example of the 
preferential attachment for scale-free networks [75] . On 
this basis, there have been attempts for the multiplex net¬ 
work evolution model [91 1561174] . 

A notable new idea for multiplex evolution model is 
that the layers coexisting in multiplex system may also co¬ 
evolve. The multiplex layers do not merely evolve together 
in time, but the evolution of layers can also become entan¬ 
gled, in the sense that evolution of one layer is dependent 
not only on the state of the current layer but also on those 
of other layers. This idea has been implemented using the 
preferential attachment scheme We illustrate the 

main idea and results following 0. 

In the coevolving multiplex network growth model, a 
node is newly added to the system in each step and con¬ 
nects to existing nodes following the growth kernel II. The 
coevolution property is implemented in the way that the 
growth kernel II a of a node’s degree in layer a depends 
not only on its degree in that layer but also its degrees in 
other layers, 

n a = f(k^M 2 \...M a \---M e) ), (5) 
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Fig. 3. (a) Schematic illustration for the coevolution model. 
Each step, a new node (orange) enters the system and estab¬ 
lishes a link in each layer. To choose the node to connect in the 
layer A, the new node refers to the nodes’ degrees not only in 
that layer A but also in the other layer B (and similarly in the 
layer B). Relative dependency to the other layer is controlled 
by the coevolution parameter e. (b) Joint degree distribution 
P(/ca, ks) for different coevolution parameter e, demonstrating 
the increasing interlayer degree correlation with coevolution, 
(c) Interlayer degree correlation coefficient p as a function of 
coevolution parameter e and shift factor a. Adapted from ref¬ 
erence 


where the superscripted layer indices are used. A simple 
example for two-layer network (Fig.[3k) is provided by the 
growth kernels in the form of coupled linear preferential 
attachment, given by [9]: 

n A oc [(1 — e)(k A + n) + e{k B + &)]> 

IIb oc [e{k A + a) + (1 - e){k B + a)], (6) 

where the constant e is the so-called coevolution param¬ 
eter and a is the constant determining the layer’s native 
degree exponent. The degree distribution of each layer is 
shown to be independent of the value of e and becomes 
p{k) ^ ^ (3+a) f or layers. However, the degree cor¬ 
relation between layers is strongly affected by the coevo¬ 


lution parameter e (Fig. [3 Jd). Specifically the interlayer 
degree correlation coefficient p defined by equation (4) is 
obtained asymptotically as 


6e + a 
6e T 2d 


(7) 


for a > 0 !9 . The interlayer degree correlation increases 
with the strength of coupling between the layers’ growth, 
quantified by the coevolution parameter e in this model 
(Fig. 3c). 

The model proposed in 156 is similar, and was later 
extended also to the case of nonlinear growth kernels m, 
showing how the structural instability might emerge in 
growing multiplex networks. 


3 Problems with cooperative layer-coupling 

In many multiplex systems, nodes can be influenced by co¬ 
operative activity from multiple layers of the system. Of¬ 
ten the overall synergistic effect is non-additive and non¬ 
linear in individual layer’s effect, leading to a situation 
that conventional single-layer network framework could 
not address accurately. It therefore has been an outstand¬ 
ing topic in multiplex network studies. In this section we 
present an overview of the class of problems regarding the 
structural and dynamical properties of multiplex networks 
in which the layers are coupled in cooperative manner. 
Problems belonging to this class include (i) the systems 
in which the proper functioning of one layer is dependent 
on the proper function of another layer and vice versa , of¬ 
ten framed in terms of interdependent networks EE and 
(ii) the processes in which a node’s state is affected by syn¬ 
ergistic influence from the states of multiple layers, such 
as multiplex threshold cascade [6], 76]. In what follows, we 
will make concise discussion focusing on how to general¬ 
ize traditional network theory to multiplex networks with 
cooperative layers and what the main effects of coopera- 
tivity in layer-coupling are, with the examples of mutual 
percolation, structural robustness, and cascading dynam¬ 
ics. 


3.1 Mutual percolation 

Being connected (or connectivity^ for short) is a minimal 
condition for the functionality ofa system. Percolation is 
a classic statistical physics problem concerning the global 
connectivity of a system. One way to generalize the con¬ 
cept of connectivity in multiplex systems with many layers 
is to require the simultaneous connectivities across differ¬ 
ent layers, leading to the notions of mutual connectivity 
and mutual percolation EHEE Motivations of such defi¬ 
nition may come from the interdependent networks of crit¬ 
ical infrastructure (37 and systems with multiple critical 

3 In graph theory, the term connectivity is used to specifically 
refer to the minimum number of nodes that need to be removed 
to disconnect the remaining nodes from each other. Here we 
are using the term in a loose sense. 
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Fig. 4. (a) Illustration for the notion of mutual percolation. On the left, even though each layer is well connected, the multiplex 
does not have a mutual component. On the right, the multiplex gains a big mutual component (marked in orange), after one more 
link is added to each layer, (b) The size of giant mutual component as a function of the layer mean degree z for multiplex ER 
networks with 1, 2, 3, 5,10 layers. The dot denotes the size of jump at the transition. (Inset) The scaling plot of M — M c vs. 
z — z c above the transition, showing the different scaling exponent /3 = 1 for i — 1 and ft = 1/2 for t — 2. (c) Mutual percolation 
of the partially-multiplexed ER networks with multiplex fraction q — 0.1 and q = 0.7, displaying continuous and discontinuous 
transitions, respectively. (Inset) The phase diagram in (g, z)-plane with mutual-percolated and unpercolated phases separated 
by lines of continuous (dotted) and discontinuous (solid) transitions, joining at the tricritical-like point (dot). 


resources m- It is a type of cooperative layer-coupling 
in the sense that each and every layer needs to operate 
concurrently for the proper functioning of the system as 
a whole. 

To formalize the problem of mutual percolation, we 
first define the mutually-connected component (or the mu¬ 
tual component for short) as the set of nodes in which each 
pair is connected within each and every layer simultane¬ 
ously [151177] (Fig. 4a). The mutual component which is 
extensive in network size, if it exists, is called the giant 
mutual component. The notion of mutual percolation was 
first introduced in [15] while describing the iterative cas¬ 
cade of failures in interdependent networks. Key finding 
of [15] was the discontinuous transition in the size of giant 
mutual component at the critical fraction of random node 
removals, colloquially referred to as the abrupt collapse 
of the system. Such a discontinuous mutual percolation 
transition is a strong departure from single-layer network 
percolations, which is a major reason for much attention 
from statistical physicists studying network theory [20 . 

Later the problem was re-casted in purely structural 
terms in m by noting that at the end of the cascade of 
failure process all the remaining (unfailed) nodes should 
form mutually-connected components. This formulation 
allows a more easily-accessible analytical approach 


1751 by extending the existing methods from single-layer 
network theory, which we will take as the basis of the pre¬ 
sentation that follows. The analytic approach put forward 
in HE which explicitly tracks the back-and-forth cascade 
steps of failure, is slightly more involved, yet both ap¬ 
proaches yield the equivalent answer when it comes to the 
giant mutual component. 

Mutual components in a given network can be found 
algorithmically by straightforward iterations towards the 
requirement of the mutual connectivity he More efficient 
algorithms for finding the mutual components have also 
been proposed [791180] . 

Analytically, the size of giant mutual component of 
multiplex network can be obtained for networks with locally- 
treelike undirected uncorrelated layers [571771178], as fol¬ 
lows. One first consider the probability w a that the node 
reached by following a randomly chosen a-layer edge does 
not belong to the giant mutual component. These proba¬ 
bilities satisfy the coupled self-consistency equations, 

n (s) 

k a b=l,bj^a 

with z a is the mean degree of a-layer and a = 1,... , A 
Terms inside the parentheses give the probability that at 
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least one neighbor in a- and all other layers, respectively, 
of the node belongs to the giant mutual component, the 
product of which becomes the probability that the node 
itself belongs to the giant mutual component under the 
condition of mutual connectivity. By averaging this prob¬ 
ability over k weighted by the factor k a P(k.)/z a , giving the 
probability that the multiplex degree of the node reached 
by following a randomly chosen a-layer edge is k, one 
can end up with equation (8). The size of giant mutual 
component is obtained by the probability M that a ran¬ 
domly chosen node belongs to the giant mutual compo¬ 
nent, which by the same reasoning is given by 

M = Y J P{ k) IT 1 -«£“)> (9) 

k a= 1 

with Wa s being the smallest solution (in the physical range 
[0,1]) of equation (8). Note that this approach reduces to 
that of the ordinary percolation for single-layer networks 
for t = 1. 

Specific application to multiplex networks with Erdos- 
Renyi (ER) layers illustrates the key features. For multi¬ 
plex ER networks with i layers of equal mean degree z, 
equations (8) and (9) reduce to a single equation for M as 

M = (l-e~ zM Y, (10) 

which undergoes a discontinuous transitions from unper¬ 
colated (M = 0) to mutual-percolated phase (M > 0) at 
the critical layer mean degree, z c (£) (Fig. 4b, main panel). 
For duplex ER networks (£ = 2), the transition point is 
z c ( 2) = 2.455407..., which is significantly higher than 
that of ordinary (single-layer) percolation, z c ( 1) = 1, and 
the jump in the giant component size at the critical point 
M c ( 2) = 0.511699... HMZZ]. 

The discontinuous mutual percolation transition shows 
some unusual properties compared to the typical first- 
order phase transition. Firstly, the transition does not ex¬ 
hibit hysteresis. Transition occurs at the same point as the 
mean degree is either increased from below or decreased 
from above the transition point, providing the basis for the 
equivalence between the two approaches of references m 
[77]. The nature of discontinuous mutual percolation tran¬ 
sition was further investigated in reference EH], showing 
that the transition is of hybrid type (or mixed order), with 
discontinuous jump in the order parameter (called the gi¬ 
ant viable cluster size in Ref. US]) followed by the critical 
scaling, 

M-M c ~(z-z c f, (11) 

above the transition point. For duplex ER networks, the 
critical exponent (3 = 1/2 (Fig. 4b, inset), as in the boot¬ 
strap and the k- core percolation on single-layer networks 
[8lj l82 ] . The average size of finite mutual components, an 
analogous quantity to the susceptibility in ordinary per¬ 
colation, does not diverge at the transition point. Instead, 
the critical scaling above the transition point is attributed 
to the diverging scale of avalanches upon the node removal 
as the transition point is approached from above m • 


3.1.1 Partial multiplexity 

Situation may occur that some of the nodes in a multi¬ 
plex do not participate in all the layers, a case that might 
be called the partial multiplexity. This generalization has 
been studied in both approaches [771183] . Let us consider 
the multiplex network in which the fraction q of nodes 
participate in both layers (multiplex nodes) but the re¬ 
maining 1 — q fraction of nodes in each layer are specific 
to that layer. The generalized mutual component is de¬ 
fined as the set of connected nodes in which every multi¬ 
plex node is connected in both layers simultaneously. In 
general the size of generalized mutual component in differ¬ 
ent layers can be different. For duplex ER networks, the 
sizes of generalized mutual components in the two layers 
are obtained by generalizing equation (10) with i = 2. It 
reads [77 

M a = (1 — e~ ZAMA )(l - qe ~ ZBMB ), 

M b = (1 -qe~ ZAMA )(l - e~ ZBMB ), (12) 

where M a and z a are the generalized mutual component 
size and the mean degree of a-layer. One can see that 
Ma 7 ^ Mb when za ^ zb in general with q ^ 1. A key 
feature of this model is that depending on q and Zi s the 
nature of transition can change from discontinuous to con¬ 
tinuous, exhibiting the tricritical-like point [77 (Fig. 4c). 

In the framework of interdependent networks, more 
general scenario was considered [7713], such that in 74- 
layer a fraction qA of nodes are dependent on a node in 
5-layer (and vice versa). The dependency relation can also 
be unidirectional, making it more general than the sim¬ 
ple multiplex networks. It can be treated analytically in 
similar way and gives qualitatively similar picture as the 
previous partial multiplex case. Through this setting with 
za = zb-> the critical behavior of the giant generalized mu¬ 
tual component was obtained, as Ma (qA - q\Y with 
(3 = 1 along the line of discontinuous transitions whereas 
(3 = 1/2 when qs = q^ B is fixed, with (q%q* B ) denoting 
the point at which the lines of discontinuous and contin¬ 
uous transitions join [83] . A phenomenological analogy to 
liquid-gas transition exhibiting the critical point was pro¬ 
posed in [83], but further understanding of the nature of 
symmetry-breaking in this problem is required to estab¬ 
lish a more complete correspondence. Generalizations of 
partially interdependent networks with both dependency 
and connectivity links have been studied in [84lf85] . 

3.1.2 Network of networks 

In the framework of interdependent layers, one can further 
extend to the general case of the network of interdepen¬ 
dent layers, dubbed as “the network of networks” (NoN) 
and studied for a number of different topology I86ll87l l88 ] . 
A notable result for NoN is that when the interdepen¬ 
dency between networks (layers) is such that there is a 
one-to-one dependency relation between all the nodes in 
each pair of connected layers, the mutual component does 
not depend on the topology of NoN [89]. It turns out that 
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the mutual component of such NoN coincides with that of 
the fully-connected NoN constructed with the same set of 
layers, z.e., the simple multiplex network. As an instance 
in which the interdependency between layers is not simple, 
the so-called configuration model of NoN was studied [90] . 
In this model the nodes are assigned a given number Q a 
of interlayer links in each layer a, which are connected in 
the configuration-type algorithm. A peculiar behavior is 
found by using the message passing calculation, that the 
NoN can undergo a series of multiple transitions when the 
interlayer degrees Q a are distributed heterogeneously m- 


3.1.3 Spatially embedded multiplexes 


Many complex systems existing in the real world are em¬ 
bedded in two dimensional (2D) Euclidean space to some 
extent [91], as in the cases of multiplex transportation sys¬ 
tem or interdependent critical infrastructure. The mutual 
percolation problem for spatially embedded multiplexes 
poses theoretical challenge as well to cope with the correla¬ 
tion imposed by the spatial constraints. First, the mutual 
percolation of duplexes with diluted 2D square lattice was 
studied [92], claiming that the mutual percolation transi¬ 
tion becomes continuous and the order parameter expo¬ 
nent 1 could be larger than that of ordinary percolation, 
suggesting that the mutual percolation transition can be 
less abrupt than ordinary percolation transition, in sharp 
contrast to the cases of random networks. Specifically in 
2D, it was estimated that (3 = 0.171(2) for duplex square 
lattices, compared to /3 S i ng i e = 5/36 = 0.1386... for a 
single-layer 2D square lattice. This claim was challenged 
in reference [93], arguing that the exponent /3 should have 
the same value in both duplex and single-layer lattices. 
Lacking the exact solution for the multiplex lattices, how¬ 
ever, the conflict has not been resolved completely as yet. 

Later, a variant of multiplex lattice model, called the 
interdependent lattice networks, in which the node in one 
lattice makes a dependency link to a node in the other 
lattice randomly chosen within a certain distance r [94] . 
For r = 0 the model reduces to the multiplex lattice model 
of [921193] . As r —)> oo, the two lattices are coupled by com¬ 
pletely random dependency links, thereby one can expect 
the transition to be discontinuous as in the multiplex ran¬ 
dom networks. In between there might exist a point r* 
separating the two behaviors. For two-layer interdepen¬ 
dent lattices it was estimated that r* ~ 8 (in lattice unit) 
[94] . For r < r*, the transition is continuous, whereas it 
becomes discontinuous for r > r*. More recently, it was 
shown that for r oo, the mutual percolation transition 
in even partially-interdependent lattices is always discon¬ 
tinuous as long as the fraction of dependent nodes q > 0, 
that is, any nonzero fraction of dependent nodes can drive 
the transition to be discontinuous [95] . This result was in¬ 
terpreted that the spatially embedded multiplex networks 
can be extremely more broadly vulnerable to abrupt (dis¬ 
continuous) collapse than random networks. 



Z 



Z total 


Fig. 5. Mutual percolation of duplex ER networks with (a) 
interlayer degree correlation of three different kinds, the MP, 
uncorrelated, and MN couplings, and (b) link overlap with dif¬ 
ferent overlap fraction r = 0.2, 0.4, 0.6. (Inset) The size of jump 
M c at the transition as a function of the overlap fraction r, 
which decreases gradually to zero as r 1. Adapted from 
references m (a) and EP (b). 


3.1.4 Effect of correlated multiplexity 

The effect of correlated multiplexity on mutual percola¬ 
tion has been addressed in terms of the interlayer degree 
correlations [35ll57U68ll96j and the link overlap [59j[60j[6T] . 
The most commonly-observed interlayer degree correla¬ 
tion would be a positive interlayer degree correlation, the 
effect of which has been the major subject of studies [35l 
|68][96]. The main result of those studies is that the positive 
interlayer degree correlation (also termed the intersimilar¬ 
ity in Ref. [35]) makes the multiplex system more robust 
to random failure p5l!57U68j , in the sense that the mutual 
percolation transition point becomes lower (Fig. 5a). The 
study has been extended to include negative interlayer de¬ 
gree correlation as well 153 • More recently, in the frame¬ 
work of interconnected networks a particular form of pos¬ 
itively correlated interlayer connection pattern was intro¬ 
duced as an optimal coupling architecture for robustness 
of the mutual connectivity to random failure [97 , offer¬ 
ing a consistent perspective with the earlier findings from 
the multiplex and interdependent networks. It was also 
proposed that the suggested optimal pattern might un¬ 
derlie the topological stability of the interconnected func¬ 
tional brain networks E3- To what extent such broad per¬ 
spective of topological stability could extend to the actual 
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functional stability of given multiplex systems remains an 
outstanding open problem. 

Another pervasive kind of correlated multiplexity is 
the link overlap. The overlapping link plays a distinctive 
role in mutual percolation, in that the nodes connected 
by overlapping links form a mutual component by them¬ 
selves. Therefore the link overlap would facilitate mutual 
percolation. In the extreme case where every link is over¬ 
lapping link, all the layers becomes identical and the mu¬ 
tual percolation problem reduces to ordinary percolation 
of the single layer, showing continuous transition. A pri¬ 
mary question of theoretical interest would then be the ex¬ 
istence of the nontrivial value of the fraction of overlapping 
links r across which the nature of transition changes from 
discontinuous (as in r = 0) to continuous (as in r = 1). 
This question was examined by two groups independently 
using different methodology BSUSD], only to reach conflict¬ 
ing conclusions. The two approaches were later scrutinized 
in EEL to show that they in general apply to different per¬ 
colation processes (see the next section for more details). 

To illustrate the effect of link overlap, let us consider 
the case of duplex ER networks of layers with equal mean 
degrees £ to tai among which the fraction r is the overlap¬ 
ping link. As the overlap cluster forms a mutual compo¬ 
nent by itself, it is useful to renormalize it as a supernode, 
which is connected by non-overlapping links to other su¬ 
pernodes (overlap clusters). To obtain the giant mutual 
component size, one needs to augment equation (10) by 
taking into account the sizes of supernodes, the distribu¬ 
tion of which is given by the component size distribution of 
the network formed by overlapping links, denoted R(m). 
The giant mutual component size M is then given by [59J 

EU 

oo 

M = tfoo + Yj R(m)(l - e~ mzM ) 2 , (13) 

m= 1 

where R 0c is the probability that a randomly chosen node 
belongs to the infinite-size overlap cluster and z is the 
mean degree for non-overlapping links, z = 2 to tai(l — 
r). Solving this, it was obtained that the transition re¬ 
mains discontinuous as long as the overlap fraction r < 1 
(Fig. 5b), with the jump M c at the transition vanishes 
gradually as r —>> 1 (Fig. 5b, inset); the tricritical point 
claimed in m is found absent. This means that despite 
its facilitating role, the link overlap per se is insufficient to 
change the nature of transition in multiplex ER networks. 
Yet the possibility of nontrivial tricritical point in layers 
of other topology remains open. Furthermore, it is left as a 
theoretical challenge to extend the analytic theory beyond 
the simple multiplex cases m ■ 

3.1.5 Related models: Viability, weak and k-core percolation 

In the mutual percolation problem, an implicit assump¬ 
tion was that the connectivity by itself is sufficient for 
functionality of the system. In some multiplex systems, 
the connectedness among nodes per se may not be suffi¬ 
cient but the connected component can be functional only 


when it is connected also to the resources essential for the 
functionality of the system. A representative example in 
this class is the modern city, which is supported by mul¬ 
tiple layers of resource supply network such as power and 
water supply networks. These vital resources are not pro¬ 
duced in every node; rather they are generated from a 
specific set of nodes (power stations and water sources, 
respectively). Taking this feature into account, a model 
of viability of multiplex system with multiple resource de¬ 
mands was proposed ma. which is defined as follows. 

Consider a network with ^-multiple layers, where each 
layer of the multiplex network corresponds to a certain re¬ 
source supply network. A given fraction p of nodes (called 
the resource nodes) generates and distributes resources 
essential to be viable. Only viable nodes can function 
properly and transmit resources further to their connected 
neighbors. Then, a node is viable only if it can reach, via 
other viable nodes, to a resource node in each and every 
layer. The viability V of the system is defined as the frac¬ 
tion of viable nodes, given the network parameters and 
distribution of resource nodes. To compute the viability 
V, one can set up equations analogous to those for the 
mutual percolation, equations (8) and (9) as [10] 

l-v a mp+(l-p)J2^^(l-v^) 

k Za 
i 

xlld-^t (14) 

6=1 

b^a 

where v a is the probability that the node reached by fol¬ 
lowing a randomly chosen a-type link (a = 1, is not 
viable, and 

i 

V = p + (l-p)J2P(Vl[(l-v k a “). (15) 

k a= 1 

In equations (14) and (15), the first term on the right hand 
side is the probability that the chosen node is a resource 
node and the second term gives the probability that it is 
not a resource node but is connected to the giant viable 
cluster. 

Applying to duplex ER networks with equal layer mean 
degree z, one has the single equation for V analogous to 
equation (10) as 

V = p + (l-p)(l-e- zV ) 2 . (16) 

Solving equation (16), one obtains two stable solutions for 
a range of parameters p and 2 ), suggesting the bistability 
in viability. Viability undergoes discontinuous jump as one 
of the two stable solutions loses its stability. The two so¬ 
lutions are shown to correspond to outcomes of different 
algorithms for identifying the viable components, called 
the cascade of activations (CA) and the cascade of deacti¬ 
vations (CD), respectively [10] (Fig. 6a). This means that 
the viability of a system depends also on the way it reached 
the current state, that is, the viability exhibits hysteresis 
(Fig. 6b). In the limit p —* 0, equations (14)-(16) reduce 
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Fig. 6. (a) Illustrations of the CA and the CD algorithms 
for viability. Resource nodes (squares) generate resources. If a 
node connects with resource nodes through each type of link, 
denoted by the solid and dashed lines, the node is viable (filled 
circles) and can transmit resources further to its neighbors. If 
not, the node is unviable (open circles). Shaded (green) circles 
denote the node whose state is to be updated (activated in CA 
and deactivated in CD) at each step, (b) Hysteresis curve of V 
with p 0 — 0.02. Starting from the well-connected high viability 
state, the systemic collapse (o) and the subsequent recovery 
(□) exhibit different curves. Dashed line indicates £ = 1 case for 
comparison, without hysteresis. Adapted from reference [un¬ 



to those for mutual percolation equations (8)-(10). In this 
sense the viability model can be considered as a general¬ 
ization of mutual percolation. In this limit, only the CD 
branch remains meaningful (the CA branch solution be¬ 
comes trivial, V = 0 for all z), which coincides with the 
solution for mutual percolation. 

In the presence of link overlap, richer picture emerges 
ED. Viability still exhibits bistability, corresponding to 
the outcomes of the CA and CD algorithms, respectively, 
but they do not come from the multiple solutions of a sin¬ 
gle equation but are obtained as the solutions of two dif¬ 
ferent equations describing the CA and CD algorithms, re¬ 
spectively. Interestingly, the two different approaches pro¬ 
posed for solving the mutual percolation problem with link 
overlap [551IHT)] describe the two algorithms, the CD algo¬ 
rithm by the method of 59 and the CA by that of [60] . 
This finding suggests that the multiplex system can in 
general respond to activation and deactivation processes 
in inherently different ways. 

This property has also been observed in similar perco¬ 
lation models called the weak bootstrap and pruning per¬ 
colations m The weak bootstrap percolation proposed 


in m is the same as the CA algorithm in the viability 
problem, whereas the weak pruning percolation is sim¬ 
ilar but not identical to the CD algorithm^] and it is 
found that the two percolation processes yield different 
results for multiplexes DU- k- core percolation was gener¬ 
alized to multiplex networks into k-core percolation non], 
in which k-core is defined as the largest subgraph with 
each node therein having at least k a edges within the sub¬ 
graph in each layer, with k = {k a }. Hybrid phase transi¬ 
tion akin to that of mutual percolation was found except 
for (&a,&b) = (1,1), for which the transition is continu¬ 
ous |100| . 


3.2 Robustness against attacks 

Network robustness has been a major topic of network 
theory m, formulated as various percolation problems. 
Compared to the usual percolation problem concerning 
random dilution (102 _. the robustness problem focuses on 
the response of the system to a variety of different ways of 
dilution, such as node removals according to their degrees, 
commonly referred to as the intentional attack [I03lll04j . 

For the class of multiplex networks with cooperative 
layer coupling, it is the giant mutual component that is the 
primary quantity for addressing the system’s robustness. 
The size of the giant mutual component M of multiplex 
networks with locally-treelike uncorrelated layers after the 
node removals can be obtained as the generalization of 
equations (8) and (9) by following reference [105] . It reads 

q 

k 

i 

x (i — u>y _i ) H(i — w b b )i 

6=1 
b^a 

i 

k a=1 

where the function </>(k) encodes the removal strategy in 
terms of the probability that the node with multiplex de¬ 
gree k is removed. For instance, for the uniform random 
removal of fraction / of nodes (as in usual site percola¬ 
tion), <p = f (constant). For the intentional attack based 
on the node’s total degree K = Y^f a =i ^ — K c ), 

where 0 denotes the Heaviside step function and K c is 
the cutoff total degree for the attack 153 • Several attack 
strategies and multiplex scenarios have been studied such 
as the targeted attack with <p(k) oc k a with adjustable 
parameter a on fully multiplexed networks (in the form 
of interdependent networks) [ 106 1 . on partial multiplexes 

4 It might be worthwhile to note that this nomenclature is 
somewhat unfortunate, in that the original bootstrap percola¬ 
tion [98] is defined as a pruning process, and the corresponding 
activating process is known as diffusion percolation [99] . 
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urn and on network of networks [108]. Localized attack 
on spatially embedded networks was also studied [ 109] . 

Network robustness property would also be strongly 
influenced by the correlated multiplexity pTIlfiTj . Against 
random failures, the positive (negative) interlayer degree 
correlation is shown to enhance (decrease) the robustness 
of mutual connectivity in multiplex random networks. For 
the case of total-degree based attack, the vulnerability de¬ 
pends not only on the correlated couplings but also on the 
initial density of the networks m • Such non-monotonic 
and protocol specific responses of multiplex networks to 
damage could complicate the design of robust coupled sys¬ 
tems mmm- as the optimal structure could strongly 
depend on the type of damage considered. 


3.3 Cascades and complex contagion 


Percolation problems are utterly concerned with the struc¬ 
tural or topological connectivity of the system. Although 
the network structure itself can be a strong determinant 
of the system’s function and robustness, there are many 
other phenomena in which the specific dynamical process 
occurring on top of the structure also confers important 
consequences [112]. In such phenomena, besides the net¬ 
work structure, each node is endowed some attribute in¬ 
dex (or state function) which is affected by the neigh¬ 
bors’ states. The particular way the nodes’ states affect 
one another defines the dynamical processes. In this sec¬ 
tion we survey a particular class of dynamical processes 
for which the cooperative layer coupling has been a crucial 
ingredient when generalized to multiplex systems, the cas¬ 
cade models. It includes the Watts-type threshold cascade 
models [6l l76llll3] . the sandpile models mm, and the 
Motter-Lai-type overload cascade models [58] . 

The threshold cascade model began as a model of how 
the behavioral adoption spreads over the social community 
m, and was later generalized and formulated further by 
mzunH]- a simple mathematical formulation by Duncan 
Watts mu has become a standard model, called accord¬ 
ingly the Watts model. Each individual (node) can in one 
of two states, active and inactive. An inactive node decides 
to activate if the fraction of neighbors who are already ac¬ 
tive exceeds the threshold. This activation process occurs 
as a cascade until no further activation can be made. It is 
known that in networks, depending on network parame¬ 
ters and threshold distribution the so-called global cascade 
might occur, by which a finite fraction of large (infinite, 
theoretically) network can be activated from a vanishingly 
small fraction of initially-active seed nodes mg. 

This model has been extended to multiplex networks 
[61176] . in which the social influences from different layers 
are integrated non-additively and nonlinearly. The ana¬ 
lytical approach applicable to locally-treelike uncorrelated 
layers proceeds in a similar but slightly different logic as 
that of mutual percolation (see Refs. [EllTH] for detailed 
arguments). The cascade size p from random-distributed 


initial seeds of fraction po can be calculated as 

k t 

p = Po + (1 - Po) E P(k) n B m 0 (?»)^( m . k ). 
k m=0 a= 1 

(19) 

where B^(q) is the shorthand notation for binomial dis¬ 
tribution, (^)g m (l — g) /c_m . The quantity { q is the 
fixed point of the coupled recursion relation, 

q n +1 - po + (i - po) 2_ j ~— 

k 

k a — 1 {kb} 

x E E B ™T(<?i a) ) 

m a = 0 {mj,}=0,6^a 

(20) 

b^a 

starting from every q^ = p 0 . F(m, k) is the response 
function encoding the activation rule, giving the probabil¬ 
ity that the node with m = {m a } active neighbors among 
its multiplex degree k will activate. This analytical ap¬ 
proach could further be generalized to multiplexes with 
correlated and modular layers by employing methods de¬ 
vised for the single-layer networks [11911120] . 

Two simple yet representative activation rules are con¬ 
sidered (Fig. 7a): (i) the OR rule, in which the node acti¬ 
vates once the threshold is met in at least one layer, and 
(ii) the AND rule, in which the node activates only after 
the threshold is met in all its social layers. Intuitively, the 
OR rule would facilitate the cascade, whereas the AND 
rule will impede it, which are confirmed by both the¬ 
ory and simulations [MTH] , A less intuitively clear picture 
emerges when the two activation rules are mixed in pop¬ 
ulation m ■ As the fraction of nodes following the AND 
rule increases, there exists a critical fraction above which 
the transition to global cascade becomes discontinuous, 
compared to the continuous transition below this fraction 
(Figs. 7b and 7c), followed by tricritical-like scaling behav¬ 
ior m- This is another instance in which the cooperative 
layer coupling, now driven by dynamics, induces discontin¬ 
uous transition in multiplex systems. Different multiplex 
generalization of Watts model has also been studied ma. 

Another prototypical model of cascades on networks is 
the model of cascading failure based on loads, put forward 
by Motter and Lai m- Each node i is assigned its own 
capacity, Ci = (1 + a)Li, setting the maximum load that 
the node can endure, where a is the tolerance parameter 
and Li is the load defined as the total number of shortest 
paths passing through that node i [122] . Initiated by the 
failure (removal) of a small set of nodes, the global redis¬ 
tribution of the shortest paths and thereby the load can 
induce the cascade of failure of overloaded nodes whose 
load becomes to exceed the capacity. Its distinctive prop¬ 
erty is that the cascade occurs in a nonlocal manner, un¬ 
like the threshold cascade models introduced above. The 
Motter-Lai-type models have been studied on multiplex 
networks with a particular focus on the effect of interlayer 
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4 Problems with complementary 
layer-coupling 

Different layers in the coupled transportation system pro¬ 
vide alternative or complementary means of traveling from 
one place to another. In such systems, malfunction of one 
single layer does not fundamentally alter the functional¬ 
ity of the whole system. In this section we consider the 
problems arising in such multiplex systems with comple¬ 
mentary layers. Complementary layers are coupled in con¬ 
nective wa}j^J and often the whole system can be regarded 
as a single super-network in which the network layers con¬ 
nected by interlayer links. Hence the problems in this class 
are admittedly much more straightforward to generalize 
from single-layer framework than the problems with, say, 
cooperative layers. This aspect might have played a role in 
prompting skepticism around the whole multiplex network 
framework. Notwithstanding, it is still useful to retain the 
multiplex framework in such problems for the sake of gen¬ 
erality and methodological clarity, as we shall argue below. 


4.1 Percolation 



Fig. 7. (a) An example multiplex configuration for illustrat¬ 
ing the two response rules. Assuming the uniform threshold for 
activation R = 1/2, green (dark) circles denote active nodes 
and the light (white, orange) ones the inactive nodes. The or¬ 
ange node will activate itself under the OR rule, but will re¬ 
main inactive under the AND rule, (b) The cascade size p with 
the initial seed fraction po = 0.001 and the uniform threshold 
R = 0.18 as a function of the layer mean degree, z, and the 
fraction of OR nodes, £. Red dot marks the cusp point at which 
the lines of continuous and discontinuous transition join. Three 
vertical planes with the guideline denote the slices for which 
the panel (c) is drawn, (c) The cascade size p, with the same 
parameter as (b), as a function of z, for three different £ — 0.2 
(green dotted), 0.5 (orange dashed), and 1.0 (red solid). For 
£ = 0.2 the transition to global cascade occur discontinuously, 
whereas it is continuous for the other two cases. Adapted from 
reference [76] . 


degree correlation and the optimal coupling patterns [58] . 
For instance, it has shown that the positive interlayer cor¬ 
relation is beneficial in mitigating the cascades compared 
with other coupling patterns, in line with the findings in 
mutual percolation [58] . 


Percolation of a system of i interconnected layers was first 
addressed systematically in 1 11 although there had been 

related predecessors like the graphs with colored edges [69j 
CD]. In this problem, one is primarily interested the exis¬ 
tence of the giant component, that is the extensive subset 
of nodes which are connected through links regardless of 
their types (layers). They developed a generating function 
method, by generalizing the well-established single-layer 
framework for the interconnected layers m- Methodolog¬ 
ically this framework could also find predecessors in the 
percolation of clustered single-layer networks [123] . Below 
we will present the framework applied specifically to sim¬ 
ple multiplexes [291157] . The general framework is not too 
different, and the interested readers are referred to El- 
Given the multiplex degree distribution, P(k), and its 
generating function, 

i 

Go(x) = £P(k)n^ (21) 

k a= 1 

with x = {x a }, the generating function for the remaining 
degree distribution after following a randomly chosen a- 

5 In this respect many problems in this class can also be 
formulated as the interconnected networks, especially in terms 
of mathematics. 

6 Comparison of the two seminal papers mum may reveal 
an interesting perspective regarding how the community per¬ 
ceives on multiplex problems of different kinds. Both works 
were, to our knowledge, first presented to wide public at the 
International Conference on Network Science (NETSCI) in 
May 2009 and uploaded in parallel onto arXiv that summer. 
Buldyrev et aV s work on mutual percolation was eventually 
published in Nature in 2010 m and is by now cited about 800 
times according to Google Scholar; Leicht and D’Souza’s work 
is still unpublished m, with some 50 citations thus far. 
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Fig. 8. Sizes of the giant component S (solid lines) and the 
giant bicomponent B (dotted lines) for the MP (red), the un¬ 
correlated (black), and the MN (blue) couplings of duplex ER 
networks. Note that the two curves for the MP coupling overlap 
completely. Adapted from reference US- 

layer edge is obtained as 

G^ a) (x) = —G 0 ( x ). (22) 

z a OX a 


and can span the whole network (*5 = 1) at a finite mean 
degree z* = 1.146193 ... per layer 0291. As for the giant 
bicomponent, in MP correlated duplexes its size B is ob¬ 
tained identical with the giant component size S', mean¬ 
ing the giant component itself becomes a bicomponent, 
whereas in general the giant bicomponent grows more 
slowly than the giant component m- 

4.2 Diffusion 

Diffusion is one of the simplest yet important dynamical 
processes, which has been studied extensively on networks 
]125U126| . A key determinant of diffusion dynamics on net¬ 
works is the spectral property of the Laplacian matrix, and 
in this sense diffusion dynamics is determined chiefly by 
the structure of networks. On this basis, there has been 
development of the mathematical formalism for studying 
diffusion processes on multiplex networks using the supra- 
Laplacian matrix EE]- Suppose a simple Gplex network 
of N nodes, in which particles diffuse with diffusion con¬ 
stant D a along each layers and D^ across different layers 
(Fig. 9a). Equations governing the time evolution of par¬ 
ticle densities x\ a ^ at each node i on layer a are given by 


The giant component size S of the multiplex of locally 
treelike uncorrelated layers can be obtained using these 
generating functions as 


S = l-G 0 (u), (23) 


dx, 


(a) 


dt 


N 

D a 


3 = 1 


A? 


(x 


( a ) _ A a ) 


) T y ^ D a t,(x, 


(*) _ » 


x) ')> 


6=1 

b^a 


(26) 


where u = {u a } is the probability that a node reached 
by following a randomly chosen a-layer edge does not be¬ 
long to the giant component, satisfying the coupled self- 
consistency equations, 

u a = G^\ u), with a = 1 , 2, ..., L (24) 

The size B of the giant bicomponent defined as the ex¬ 
tensive subset of nodes connected by at least two disjoint 
paths H23I is also calculated for multiplex networks as m 


with wff being the link-weight matrix of layer a, which 
can be recasted in matrix form as [7 


dx 

dt 


= Ex. 


(27) 


with x = (xp ,..., x $,..., xp *,..., x $) T and C being 
the supra-Laplacian matrix of the multiplex. For duplexes, 
the supra-Laplacian matrix can be written explicitly as 


B = 1 — £? 0 (u) — ^(1 — u^ZaG^tu). (25) 

a 


( d a l a + d ab i -d ab i \ 

\ —DabI D b L b + D ab I J ’ 


(28) 


The last term gives the difference between the giant com¬ 
ponent S and bicomponent B. The condition for existence 
of the giant component S > 0 and that for the giant bi¬ 
component B > 0 coincide and is given as that the largest 
eigenvalue of the Jacobian matrix J of equation (24) at the 
trivial solution u = {1,..., 1} be larger than unity [57] . 

Applying to duplex ER networks with layers of equal 
mean degree z, a number of peculiar behaviors of cor¬ 
related multiplex random networks have been observed 
(Fig. 8). First, the MP interlayer degree correlation was 
shown to promote the percolation to the extreme, in that 
the percolation threshold decreases all the way to z c = 0. 
On the other hand, in the MN correlated duplexes, the 
percolation is significantly delayed as z c = 0.838587..., 
compared to z c = 1/2 for the uncorrelated coupling, but 
once it occurs the giant component grows more quickly 


where L A and L B are the Laplacian matrix of each layer, 
defined by L a = S a — W a , where W a is the weights matrix 
of layer a and S a is a diagonal matrix with the elements 
(. S a )ii = [17!. Of particular importance is the 

smallest nonzero eigenvalue A 2 of the (supra) Laplacian 
matrix C which sets the diffusion timescale r as r = I/A 2 , 
which characterizes how fast the system could relax to the 
stationary state (corresponding to Ai =0). 

By investigating the eigenvalue spectrum of £, equa¬ 
tion (28), a number of general conclusions have been drawn 
[7], some of which we highlight below. First, the diffu¬ 
sion timescale r or equivalently A 2 undergoes a qualitative 
change at certain threshold interlayer diffusion strength 
Dab (Fig. 9b): if D AB is below the threshold, the dif¬ 
fusion timescale is completely governed by the interlayer 
diffusion, as r = 1/(2 D AB ). On the other hand, for D AB 
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Fig. 9. (a) Schematic illustration for the diffusion process on a 
duplex network, with the intralayer diffusion constant, Da and 
Db , and the interlayer diffusion constant, Dab ■ (b) Compari¬ 
son between the second smallest eigenvalues A 2 of the different 
Laplacians, as a function of the interlayer diffusion constant 
Dab (denoted as D x here) with Da — Db — 1* Note that 
the diffusion timescale r is given by r — 1/A 2 . Panel (b) is re¬ 
produced from reference [7] with permission from the authors 
and the publisher. Copyright (2013) by the American Physical 
Society. 


above the threshold, the timescale becomes dependent on 
the details of multiplex coupling, but it is always that the 
diffusion on the multiplex is faster than the diffusion in 
the slowest layer of the two. Only for sufficiently large 
Dab 1 5 the diffusion on multiplex can become faster 
than the diffusion on any individual layer in isolation, re¬ 
ferred to as the super diffusive property in [7pj Further 
investigation of the supra-Laplacian matrices of intercon¬ 
nected networks has identified the structural transition 
underlying the non-analyticity in diffusion timescale [16], 
which was later shown to be related with the reducibility 
transition of the supra-Laplacian matrices Ml- 

7 It is another misleading terminology, in that the term su¬ 
perdiffusion originally refers to the superlinear scaling of mean- 
square displacement in time in diffusive processes. Here the 
superdiffusive property denotes merely an enhancement of ef¬ 
fective diffusion constant, rather than a change in dynamic ex¬ 
ponent. It is an intriguing open question whether scaling prop¬ 
erties of diffusion dynamics and random walks such as spectral 
dimension mi could be modified by multiplex coupling. 


Question of whether or not would it be possible to re¬ 
duce a multiplex network into an equivalent single-layer 
network by layer aggregation was addressed also in terms 
of the supra-Laplacian matrix spectrum [128] . in which 
it was claimed that for many cases the multiple layers 
of a multiplex could be reducible without significant loss 
of information. This result, however, may originate from 
the linear nature inherent to diffusion dynamics that the 
supra-Laplacian describes, and its applicability to general 
nonlinear dynamical problems is yet to be scrutinized. 
Meanwhile, as the Laplacian spectrum is closely related 
also with the stability of synchronized state [12911130] . it 
was used in finding optimal coupling parameter for syn- 
chronizability in multiplex, interconnected systems ED- 
Several other diffusion-related problems have been stud¬ 
ied in the literature, such as the reaction-diffusion pro¬ 
cesses msm- Readers interested in further details are 
referred to a recent review on this subject 133 . 

Non-diffusive transport such as the ballistic transport 
along the shortest paths can be relevant in such problems 
as the network routing and congestion control based on 
traffic load (or betweenness) [12211134] . A pioneering work 
pointing out the relevance of multiple layers in such prob¬ 
lems m suggested the framework of what was called 
the layered network of the logical layer defined by the 
routing table and the physical layer over which the real 
transport occurs, and applied it to real traffic time table 
data of public transportation system to simulate real traf¬ 
fic patterns m- More recently, the multiplex framework 
has been applied in a more conceptually straightforward 
manner to the study of interconnected transportation sys¬ 
tems m- Metropolitan transportation network is mod¬ 
eled as a coupled (interconnected) network of spatially- 
embedded layers. Depending on the origin-destination dis¬ 
tribution, the optimal coupling strength may exist, for 
which the trade-off between time efficiency and conges¬ 
tion level can yield maximum utility of the system m- 
The navigability of multiplex transportation system upon 
random failures was studied [35 under both diffusive and 
shortest-path based transport scenarios. 


4.3 Epidemic spreading 

Epidemic spreading is another classic topic of dynami¬ 
cal processes on networks with its own extensive litera¬ 
ture [137] , It models how the infectious entities, be it a 
disease or a fad, spread over the network. The way the in¬ 
fection is transmitted to other nodes over the links defines 
a specific model. As the infection proceeds, it is not un¬ 
common that agents (say individuals) in the network can 
be exposed to infectious entities (say human pathogens) 
through many different channels of interaction (say dif¬ 
ferent social network layers like family or co-worker ship) 
which can have vastly different spatiotemporal character¬ 
istics. To handle such scenario, the multiplex framework 
would not only be more straightforward conceptually but 
also more transparent methodologically than the single¬ 
layer counterpart, which in itself should be merely effec¬ 
tive and approximate. 
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Works on classical susceptible-infected-removed (SIR) 
and susceptible-infected-susceptible (SIS) epidemic spread¬ 
ing models on interconnected networks fT38lll39H140lH4T] 
should be noted as the predecessors of bona fide multi¬ 
plex models. These models on interconnected networks 
could be studied by direct generalization of the single-layer 
methods, similarly to percolation and diffusion problems. 

Contagion process over multiplex social network with 
multiple channels (layers) of contagion as well as inter¬ 
layer contagion following SIS-type dynamics was formu¬ 
lated by means of the contact-based discrete time Markov 
chain [52]. It has shown that the epidemic threshold in 
such systems is completely governed by the layer with 
the largest maximum eigenvalue of the contact probability 
matrix and that simply aggregating different layers into a 
single network cannot describe the process accurately |5j. 
Analogous conclusion was drawn from the studies on sim¬ 
ilar models [142U143] . 

Often the coupled multiple layers contributing to spread¬ 
ing process are subject to some level of spatiotemporal 
separation. Consider the online and offline communication 
channels for information spreading, for instance. The con¬ 
cept of layer-crossing overhead (or layer-switching cost) 
was introduced as a new dynamic ingredient to account 
for this aspect in information spreading process over mul¬ 
tiplex social networks [ 144] . It implies the difference in 
the rates of across-the-layer and along-the-layer infections, 
leading to path-dependent transmissibility even for the 
same link (Fig. 10a). Studying an SIR-type model with 
such path-dependent transmissibility |144| . it was found 
that the layer-crossing overhead confers nontrivial effects 
on spreading dynamics, viz., the epidemic threshold be¬ 
comes dependent on the layer-crossing overhead in non¬ 
monotonic ways (Fig. 10b). Novel analytical method tai¬ 
lored to deal with the path-dependent transmissibility had 
also to be developed |144| . In this sense, it can be said that 
even the simple contagion model (like classical SIR model) 
can turn out not-so-simple on multiplex networks. 

Another situation for which multiplex framework is 
useful is the case of two or more infectious entities spread¬ 
ing through its own network layer over a multiplex pop¬ 
ulation. Particularly interesting is when the multiple in¬ 
fectious entities interact, either compete or cooperate, for 
spreading. The cases of competitive multiplex spreading 
were studied more actively, in the context of interplay be¬ 
tween awareness and disease spreading [145U1461114711148] , 
interaction between two diseases coexisting in a host popu¬ 
lation such as AIDS and tuberculosis [149]. and spreading 
of competing memes m, providing the conditions for 
coexistence or exclusive outbreak of different diseases in 
a population. The coinfection model is worthwhile 

to note even though it is studied on single-layer networks, 
as it provides a conceptual framework on how the cooper¬ 
atively coupled epidemics (syndemics) could be modeled 
in multiplexes. It was found for the coinfection model the 
epidemic transitions can occur in a discontinuous manner 
together with hysteresis hsb, consonant with the find¬ 
ings from mutual percolation and viability problems in 
Section 3. 



6X 


Fig. 10. (a) Schematic illustration of the SIR-type model with 
layer switching cost on a duplex network. Transmissibility T for 
the infection along a link is determined by the types of both 
incoming (denoted by second subscript) and outgoing (first 
subscript) infection channels (layers), (b) Color-coded plot of 
the epidemic threshold A c as a function of the layer density 
disparity parameter Sz , defined by za,b — (l±£z)zo/2, and the 
infection rate disparity parameter S A, defined by Amtra,inter = 
(1±£A)A, on duplex ER networks with zo = 2.5. The horizontal 
dashed line (black) is the boundary above which the epidemic 
threshold decreases monotonically with 5 A. Green (blue) line 
indicates the loci with the lowest (highest) value of A c for given 
5z , displaying non-analyticity. Adapted from reference pm- 

5 Other types of layer coupling 

5.1 Competitive coupling, antagonistic layers 

Another important type of layer couplings in multiplex 
systems is the competitive coupling of antagonistically in¬ 
teracting layers [15211153] . It can also be regarded a sub- 
type of cooperative coupling in that the layers are inter¬ 
locking one another functionally, but we shall treat it sep¬ 
arately to emphasize its unique features. 

The concept of mutual percolation has been twisted 
into the percolation of antagonistic layers [15211154] (which 
might well have been called the exclusive percolatiorQ . In 
this model, for two-layer multiplexes, the percolating clus¬ 
ter of each layer is defined in the following way: a node i 

8 Not to be confused with the explosive percolation \m ■ 
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belongs to the percolating cluster of layer A if it has at 
least one neighbor in layer A belonging to the percolating 
cluster of layer A, but at the same time, has no neighbors 
in layer B belonging to the percolating cluster of layer B. 
{vice versa for layer B.) The equations for the size of per¬ 
colating clusters in each layer S a can be set up by following 
similar logic as equations (8) and (9), as follows. First, the 
probability w a that the node reached by a randomly cho¬ 
sen a-layer link does not belong to the percolating cluster 
of layer a satisfy the coupled self-consistency equations, 


w { 




w 


ka lj w k b 


1 -G[ a \w a ) G ( 0 °’(w b ), 


(b), 


(29) 


for distinct a, b G {1,2}, and the last equality holds when 
P{ k) factorizes (that is, the degrees of a node in differ¬ 
ent layers are uncorrelated). (x) and G < f J \x) are the 
generating functions defined as in equations (21) and (22) 
but for each layer a. Sizes of percolating clusters S a are 
given by [ 152] 




E P (k)(l-«a“K l 

k 




Go\w b ), 


(30) 


with the same convention as equation (29). Solving equa¬ 
tions (29) and (30) for a number of different layer topology, 
a rich phase diagram was obtained. Notably the bistable 
phase accompanying discontinuity and hysteresis in sta¬ 
ble solutions was found, which occurs when the densi¬ 
ties of both layers become sufficiently large |152j and sus¬ 
tains even with only a small fraction of antagonistic nodes 
|154| . These results might be compared with those found 
in epidemic spreading with competitive layers pT45iri46l 
114911150] for their similarity and differences (see Sect. 4.3). 
It is worth to remark that it has not yet been explicitly 
demonstrated how the different percolating clusters cor¬ 
responding to solutions of equations (29) and (30) in the 
bistable phase could materialize from specific percolation 
processes, which is related to the question of the topo¬ 
logical and algorithmic definition of the giant exclusively- 
connected component and its uniqueness. Finally, a 
conceptually-related study on two interconnected networks 
competing for spectral dominance [153] is worth a men¬ 
tion, in which successful interconnection strategy was in¬ 
vestigated for a number of cases of different layer topology. 
Overall, compared to the previous two classes, much lit¬ 
tle attention has been paid to the competitively coupled 
layer problems thus far, remaining a fertile ground for fu¬ 
ture works. 


5.2 Directional coupling, hierarchical layers 

Layers in a multiplex system can also have directional or 
hierarchical, rather than bidirectional or mutual, influence 
from one layer to another. Although this type of coupling 


has yet received little attention in the multiplex network 
literature, there is another active research field, the tempo¬ 
ral network, in which the concept of directionally coupled 
layers is imposed naturally by the causality constraint due 
to time ordering. It is therefore expected that theoretical 
framework and understanding therefrom can be applicable 
to this type of multiplex problems [4]. The temporal net¬ 
work is itself an active branch of contemporary network 
theory, with growing body of literature. Readers interested 
in more details are referred to the recent comprehensive 
review and compendium on the subject [2X1(22] . 


6 Conclusion and outlook 

In this Colloquium paper, we have presented an overview 
of recent development in multiplex networks from the view¬ 
point of statistical physics. In doing so, we have put con¬ 
sistent emphasis on the importance of bringing the proper 
context of multiplexity into the problem and its solution. 
Examples of multiplex systems and multiplex structural 
measures are introduced broadly yet concisely. Various 
statistical physics problems on the structure and dynam¬ 
ics of multiplex networks are categorized by the functional 
context of layer coupling. Two major categories, the prob¬ 
lems with cooperative layers and those with complemen¬ 
tary layers, are discussed in detail and highlighted for their 
distinction. 

No review papers can be completely exhaustive, and 
this paper is no exception. In a due manner, we have cho¬ 
sen deliberately to do a highly focused review centered on 
the strict definition of multiplex networks and the prob¬ 
lems thereon. Closely-related variety of multilayered net¬ 
works are drawn into in-depth discussion only when it is 
largely indispensable. Most of concepts and methodology 
can readily be applicable to related multilayered networks 
with minimal customization. Readers are referred to [uca 
[3j to fill the possible gap, however. We have also limited 
the discussion for the most part on multiplex random net¬ 
works for the sake of simplicity. Multiplexes of more het¬ 
erogeneous layers such as scale-free networks could furnish 
a far richer repertoire of new behaviors due to the interplay 
between the structural singularity and multiplexity hmzsi 
1561, a fertile ground for further study. Last but not least, 
several other dynamical processes have also been studied 
in multiplex framework, detailed discussion of which was 
regretfully omitted in this paper. Notable examples in¬ 
clude the Boolean dynamics on multilevel networks m, 
synchronization in interconnected networks 
U evolutionary game dynamics [16211163] . and opinion 
dynamics on multiplex coevolutionary networks [164] . 

Despite some skepticism over the multiplex framework 
in favor of simpler and more-established single-layer one, 
the compelling viewpoint of this paper is that the multi¬ 
plex framework is not only useful but crucial in elevating 
the level of our understanding of complex systems. Novel 
phenomena unforeseen in traditional single-layer frame¬ 
work can arise as a consequence of the coupling of network 
layers, especially when the layers are coupled in coopera¬ 
tive or competitive manner, in which cases the multiplex 
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is not reducible into an equivalent single-layer system. 
It only remains to be seen what other new physics can 
emerge in this framework. In this prospect, useful lessons 
might be learned also from the vast literature on emergent 
phenomena in coupled statistical mechanical models m 
I166ill67j . as well as in layered physical systems such as 
layered cuprates in high-T c superconductors [ 168] , layered 
materials in multiferroics [T69] . and layered graphenes m 
We would like to close this Colloquium paper with an 
idiosyncratic list of a few research questions warranting 
immediate attention, (i) Proper context-dependent defini¬ 
tion and analysis of network measures such as the shortest 
path and betweenness centrality need to be established; 
(ii) There still remain major network concepts anticipat¬ 
ing to be applied and generalized to multiplexes. A prime 
example is the controllability |1711I1721I173U174] : (iii) No 
understanding in statistical physics is complete without 
systematic cataloguing and classification of possible uni¬ 
versality and associated critical phenomena. Compared to 
single-network counterpart [175] . the current level of our 
understanding is far from satisfactory. Identification of the 
minimal couplings (in the renormalization group sense) 
relevant to the characteristic discontinuous transitions in 
multiplex systems would be a seminal stepping stone in 
this regard. 
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